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10 Abstract. Developing accurate and tractable mathematical models for 

11 partially molten systems is critical for understanding the dynamics of mag- 

12 matic plate boundaries as well as the geochemical evolution of the planet. 

13 Because these systems include interacting fluid and solid phases, develop- 

« ing such models can be challenging. The composite material of melt and solid 

is may have emergent properties, such as permeability and compressibility, that 

i6 are absent in the phase alone. Previous work by several authors have used 

I? multiphase flow theory to derive macroscopic equations based on conserva- 

18 tion principals and assumptions about interphase forces and interactions. Here 

19 we present a complementary approach using homogenization, a multiple scale 

20 theory. Our point of departure is a model of the microstructure, assumed to 

21 possess an arbitrary, but periodic, microscopic geometry of interpenetrat- 

22 ing melt and matrix. At this scale, incompressible Stokes flow is assumed to 

23 govern both phases, with appropriate interface conditions. 

24 Homogenization systematically leads to macroscopic equations for the melt 

25 and matrix velocities, as well as the bulk parameters, permeability and bulk 

26 viscosity, without requiring ad-hoc closures for interphase forces. We show 

27 that homogenization can lead to a range of macroscopic models depending 

28 on the relative contrast in melt and solid properties such as viscosity or ve- 

29 locity. In particular, we identify a regime that is in good agreement with pre- 

30 vious formulations, without including their attendant assumptions. Thus this 

31 work serves as independent verification of these models. In addition, homog- 
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32 enization provides a consistent machinery for computing consistent macro- 

33 scopic constitutive relations such as permeability and bulk viscosity that are 

3 4 consistent with a given microstructure. These relations are explored numer- 

35 ically in a companion paper. 
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1. Introduction 

36 Developing quantitative models of partially molten regions in the Earth is critical for 

37 understanding the dynamics of magmatic plate boundaries such as mid-ocean ridges and 

38 subduction zones, as well as for providing a better integration of geochemistry and geo- 

39 physics. Beginning in the mid 1980's there have been multiple derivations of macro- 
« scopic equations for magma dynamics that describe the flow of a low-viscosity fluid in 

41 a viscously deformable, permeable solid matrix [McKenzie, 1984; Scott and Stevenson, 

42 1984, 1986; Fowler, 1985, 1989; Spiegelman, 1993; Stevenson and Scott, 1991; Bercovici 

43 et al, 2001a, b; Ricard et al, 2001; Hier-Majumder et al, 2006; Bercovici and Ricard, 

44 2005, 2003; Ricard and Bercovici, 2003; Ricard, 2007]. The details and specific processes 

45 included, vary slightly among these model systems but all are derived using the methods 

46 of multiphase flow [e.g. Drew and Segel, 1971; Drew, 1971, 1983]. 

47 Multiphase flow techniques are well formulated in many texts, including Drew and 

48 Passman [1999]; Bremen [2005]. Typically, the two-phase medium is examined at a 

49 macroscopic scale, much larger than the pore or grain scale, and one attempts to develop 
so effective media equations based on conservation of mass, momentum, and energy. This 
si approach is reasonably straightforward and has proven useful in applications, notably di- 

52 lute disperse flows. However, they have two fundamental sources of uncertainty. First, 

53 an appropriate "interphase force" must be posited. There is often a tremendous range 

54 of mathematically valid choices for this force with little to constrain it beyond physical 

55 intuition and experimental validation. Second, the macroscopic equations derived for the 

56 partial melt problem include critical constitutive relations, such as permeability, bulk vis- 
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57 cosity or effective shear viscosity. These should depend on the microscopic distribution of 

58 melt and matrix, information that is often lost in the multiphase flow approach. Multi- 

59 phase flow does not naturally determine these relationships. As with the interphase force, 
eo these closures require estimates from scalings, numerical simulations, and experiments. 

61 In this and a companion paper [Simpson et al, 2008b], we present a complementary 

62 method for deriving effective macroscopic equations using the methods of homogenization 

63 theory. Rather than assume macroscopic equations and then seek closures for constitutive 

64 relations, we assume microscopic equations and derive the macroscopic equations. This 
es is done by a multiple scale expansion, which encodes both fine and coarse length scales 
66 into the field variables. As in all multiple scale methods, the equations are matched at 
e? each order of the small parameter and solved successively. For a useful introduction to 
es homogenization with applications see Hornung [1997]; Torquato [2002]. More rigorous 

69 mathematical treatments are presented in Bensoussan et al. [1978]; Sanchez-Palencia 

70 [1980]; Cioranescu and Donato [1999]; Chechkin et al. [2007]; Pavliotis and Stuart [2008]. 
?! Homogenization has been used extensively for flow in rigid and elastic porous media, but 

72 we believe this is the first application to the magma dynamics problem which permits 

73 viscous deformation of the matrix. 

74 This strategy has several advantages with respect to multiphase flow methods. In partic- 

75 ular, there is no under-constrained interphase force, as these effects are described precisely 

76 by boundary conditions between the phases at the micro scale. Second, and perhaps more 

77 importantly, it provides a mechanism for computing consistent macroscopic constitutive 

78 relations for a given microscopic geometry. For the magma dynamics problem, it yields a 

79 collection of auxiliary "cell problems" whose solutions determine the bulk viscosity, shear 
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so viscosity, and permeability of the medium consistent with the micro-structure. We em- 

81 phasize that these macroscopic effective quantities are not volume averages of microscopic 

82 quantities. Indeed, permeability and bulk viscosity are undefined at the grain scale, but 

83 they appear as macroscopic properties through homogenization. More generally, homog- 

84 enization can extract tensor valued permeabilities and shear viscosities for anisotropic 
as media. The methods presented are also adaptable to other fine scale rheologies and 
se physics. 

a? In this work, we specifically consider the simplest case of homogenization of the mo- 
ss mentum equations for two coupled Stokes problems involving a high viscosity phase (the 

89 solid matrix) and a low viscosity fluid. This work is adapted from studies of sintering and 

90 partially molten metal alloys, [Auriault et al, 1992; Geindreau and Auriault, 1999], which 

91 in turn is based on earlier work in poro-elastic media [e.g., Auriault, 1987, 1991a; Auriault 

92 and Boutin, 1992; Mei and Auriault, 1989; Mei et al, 1996; Auriault and Royer, 2002]. 

93 For clarity, we only consider linear viscous behavior for the solid, as this may be appro- 

9 4 priate for the diffusion creep regime [e.g., Hirth and Kohlstedt, 1995a]. This assumption 

95 considerably simplifies the analysis. Extensions to power-law materials are discussed in 

96 Geindreau and Auriault [1999]. 

97 We demonstrate that, depending on the choice of scalings, we can derive homogenized 

98 macroscopic equations for three different regimes, and identify a particular regime that is 

99 consistent with existing and commonly used formulations such as McKenzie [1984] and 

100 Bercovici and Ricard [2003] . This provides independent validation of these other systems. 

101 We also discuss the strengths and weaknesses of homogenization and identifies some open 

102 questions. We recognize that the derivation is somewhat technical but have attempted to 



DRAFT 



September 4, 2009, 4:07pm 



DRAFT 



SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS X - 7 

103 make the overall approach as accessible as possible with the hope that other researchers 

104 will extend these methods to related problems (e.g., including surface energies and more 
us complex rheologies). 

we The second paper is more practical and provides specific computation of cell-problems 

10? to calculate consistent permeabilities, bulk-viscosities and effective shear viscosities for 

108 several simplified pore geometries. In particular we provide a derivation for the bulk 

109 viscosity and demonstrate that, for a purely mechanical coupling of phases, it should scale 
no inversely with the porosity; a relationship we conjecture is insensitive to the specifics of 
in the microscopic geometry. Such an inverse relationship has been suggested before [e.g. 
n2 Schmeling, 2000; Bercovici and Ricard, 2003]; however, this is the first rigorous derivation 
n3 from the microscale. Further implications of these results are discussed in the second 
iu paper. 

us Simpson et al. [2008a, b] are based on the PhD. thesis of G. Simpson, Simpson [2008]. 
2. Problem Description 

2.1. Macroscopic and Microscopic Domains 

He To begin the upscaling procedure, we must describe the spatial domains occupied by 

n7 each phase. We denote with symbol Q the total macroscopic region of interest, containing 

us both the melt and the matrix, with a characteristic length scale L which might be an 

n9 observed macroscopic characteristic wavelength (e.g. 1 m-10 km). Initially, we assume 

120 that within Q the matrix has a periodic microstructure. A two-dimensional analog is 

121 pictured in Figure 1. Q, the bounded gray region, is tiled with a fluid filled pore network 

122 of period £. £ is a representative measure of length scale of the grains or pore distribution, 



DRAFT 



September 4, 2009, 4:07pm 



DRAFT 



X - 8 SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

123 such as a statistical moment of the grain size distribution, and is much smaller than L. 

124 Notation for the domains is given in Table 1. 

We form the first important dimensionless parameter, e, 




e will play two important roles in what follows. First, all other dimensionless numbers 
and parameters will be expressed in powers of e. Second, we will expand the dependent 
variables in powers in e as in 

$ = $(°) + e$« + e 2 $( 2 ) + . . . (2) 

125 Of course, real partially molten rocks are not a periodic medium. Pore structures similar 

126 to those expected in peridotite appear in Figure 3. Since it is crystalline, there is some 

127 regularity, but it is closer to a random medium. While only the periodic case is treated in 
us this work, the random one is of interest and is also amenable to homogenization, Torquato 

1 2 9 [2002]. 

We divide our domain Q from Figure 1 into three subregions: 

The fluid portion of f2. 
il s — The solid portion of Q. 
T— The interface between fluid and solid in Q. 

130 We shall write equations for the melt in Qf, equations for the matrix in fl s , and boundary 

131 conditions between the two along T. 

132 We now introduce the notion of a cell. The cell, appearing in Figure 2 and denoted with 

133 the symbol Y, is a scaled, dimensionless, copy of the periodic microstructure of Figure 1. 

13 4 This is divided into a fluid region, Yf, a solid region, Y s , and an interface, 7. A simple 
DRAFT September 4, 2009, 4:07pm DRAFT 



SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS X - 9 

us three-dimensional example of such a cell appears is displayed in Figure 4. The cell should 

we be interpreted as a scaled representative elementary volume of the grain scale. It may be 

is? a single grain or a small ensemble of grains. 

us Although the connectedness of both phases is an important property, the particular 

139 microstructure of Y does not play a significant role in the form of macroscopic equations. 

140 Cell geometry does determine the magnitudes and forms of the constitutive relations 

141 appearing in the equations. This is discussed in the companion paper. 

2.2. Grain Scale Equations 

At the microscale, we assume both phases are incompressible, linearly viscous, isotropic 
fluids. The rheology for each phase is: 



a = -pi + 2/ie (v) (3) 
where e(v) is the strain-rate tensor, 



e(v) = \ (Vv + (Vvf ) (4) 



Components may be accessed by index notation: 

1 / dvi dvj \ 

<Jii = -p$ij + (v) 

142 The variable X appearing in these expressions denotes the dimensional spatial vari- 

143 able. The stress in il s for the solid phase is a s , with pressure p s and velocity v s . Similarly, 

144 the fluid has stress a* , pressure pf, and velocity in Qf. The notation for the fields is 

145 given in Table 2. 

At the pore scale, the Reynolds number is small; using the the values of Table 3, the 

value in the melt is < 0(1CT 5 ) and as low as 0(1CT 30 ) in the matrix. Therefore, we will 
DRAFT September 4, 2009, 4:07pm DRAFT 



X - 10 SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

omit inertial terms in the conservation of momentum equations. Each phase satisfies the 
Stokes equations at the grain scale; the divergence of the stress of each phase balances 
the body forces. As we are interested in buoyancy driven flow, the forces g s = —p s ge 3 
and = —pjge 3 are included. The equations are: 

V • a f + g f = in tt f V • a s + g s = in tt s (5a) 

V-v^O infij V-v s = infl s (5b) 

Conditions at the interface between fluid and solid, T, are still needed. As both are 
viscous, we posit continuity of velocity and normal stress: 

v s = v 7 , on T (6a) 
a s ■ n = af ■ n, on T (6b) 

M 6 A Boussinesq approximation has been made by taking the velocities to be continuous 

«7 as opposed to the momenta. These equations are exact in the sense that, subject to 

us boundary conditions on the exterior of Q, solving them would provide a full description 

H9 of the behavior of the two-phase medium (although it would be impractical to solve such 

wo a system at the macroscopic scale of interest). 

2.3. Scalings 

Our effective equations emerge from multiple scale expansions of the dependent vari- 
ables. The dimensional spatial variable X specifies a position within either Q s or Qf. We 
introduce two dimensionless spatial scales, y, the "fast" spatial scale, and x, the "slow" 
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spatial scale. These relate to X, and one another, as: 

y = X/* (7a) 
x = X/L = ey (7b) 

The expansion in (2) is now made more precise. All variables are assumed, initially, to 
have both fast and slow scale dependence: 

$( y ) = $(°)(ey,y) + e$«(ey,y) + eV\ey,y) + ... (8) 

151 Such an expansion captures grain scale detail in the second argument, while permitting 

152 slow, macroscopic variations in the first argument. As we take our domain to be peri- 
ls odically tiled with scaled copies of the cell Y, we assume $^^(x, y) is y-periodic at all 

15 4 orders of j. We seek equations that are only functions of x; these will be the effective 

155 macroscopic equations. 

Before making series expansions in powers of e, the equations must be scaled appro- 
priately. In addition to e, there are several other important dimensionless numbers. As 
motivation, let P s , P* , V s , V- f be characteristic pressures and velocities for the solid and 
fluid phases. We write: 

pf = pfpf p s = p s p s (9) 

v / = yfjrf v S = \/ S V S (10) 

Tildes reflect that the variables are now dimensionless and 0(1). Using these definitions, 
we non-dimensionalize (5a). We are free to scale the equations to either the slow or fast 
length scale. In this work, we scale to the £ length, though this does not affect the results. 
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On the length scale £, the force equations are: 



-p f I 



+ 
+ 



pf 

p s g£~ 



g/ = 

g s = o 



(ii) 



where e y denotes the strain-rate tensor with velocity gradients taken with respect to the 
fast length scale y. This motivates defining four more dimensionless numbers: 

p s V s 



^ ~ Pf£ 
f _ Pf9* 



Q 



TV, = 



Pf 



nt = 



P s £ 

_ p s g£ 



ps 



(12) 
(13) 



The Q's measure the relative magnitudes of the viscous forces and the pressure gradients, 
while the 7£'s measure the relative magnitudes of the body forces and the pressure gradi- 
ents. The /i's and g's remain in the equations as 0(1) constants. Three other important 
parameters are the ratio of the viscosities of the two phases, the ratio of the velocities of 
the two phases, and the ratio of the pressures of the phases: 



Ps 



V 



V = 



V s 

ps 



(14) 
(15) 
(16) 



156 A full list of dimensionless numbers is given in Table 4. 

Starting with e, V, and A4, we estimate these parameters with the data in Table 3: 



M = 0(1Q- 21 - 10 



-14\ 



V = O(K) 1 - 10 3 ) 
e = 0(l(n 7 - 1(T 2 ) 



(17) 
(18) 
(19) 



DRAFT 



September 4, 2009, 4:07pm 



DRAFT 



SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS X - 13 

For these values, M. <C e <C V. Since we expand the equations in powers of e, we relate 
all quantities to e. A quantity Q is said to be 0(e p ) if 

e p+1 < Q < e"" 1 (20) 

In terms of e, .M and V are approximately: 

M = 0(e n - e 2 ) (21) 
V = 0(e° - e- 2 ) (22) 

We emphasize that this power of e scale is less precise than a power of 10 scale. For 
example, V might be O^O 1 ), but if e = O(10~ 4 ), we would say V = 0(e° = 1) since 
V e -1 . Indeed, in one of the scaling regimes we consider, V = 0(e° = 1). 

To estimate the other parameters, we need estimates of the characteristic pressures. To 
do this, we first consider the forces on the matrix. At the macroscopic scale, the melt is 
0(1%) of the medium's volume. We thus argue that on this scale, the matrix is "close" to 
satisfying the Stokes equations; the pressure gradient, viscous forces, and gravity balance 
one another. On the large length scale L, the dimensionless form of (5a) is 



fH|g|2^(v s ) 



+ P ~^S S = 0. (23) 



Similar to Equation (11), we define 



Si - ^ (24) 

n - ^ (25) 
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For the terms to be in balance, Q S L = 0(1) and 1Z S L = 0(1). Using (24-25) and the 
definition of e, the fast length scale parameters are: 



C2 = O (e- 1 ) 
n\ = O (e 1 ) 



(26) 
(27) 



In the absence of direct pressure measurements, we assume the pressures are the same 
order, 

V = 0(1 = e°). (28) 

An argument for this is given in Drew [1983]. Briefly, since the velocities of interest are 
far less than the speed of sound, it would be difficult to support large pressure gradients 
across the phases without surface tension, a mechanism we do not include. We make an 
additional assumption that there are 0(1) non- hydrostatic pressures in both phases; if 

V = Phydro + Pnon-hydro then 

0(1 = e°). (29) 



V 



Pnon-hydro 

In the fluid, since p//p s = 0(1), a consequence of V — O(l) is 



i>f!i ! _ PsgZpfP 3 -p J'/.p n < ' 



pf 



(30) 



The fluid's force ratio is 



Q{ = ^ = [VMVQD = O {e-'MV) 



(31) 



and 



Q{ =MV 



(32) 
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Therefore, 



MV = 0(e 



10 




(33) 



Q{ = 0(e 9 - e 1 



(34) 



wo The choice of dimensionless parameters will lead to different expansions and effective 

161 equations. In the terminology of Auriault [1991a, b]; Geindreau and Auriault [1999]; 

162 Auriault et al. [1992], we can derive one of several outcomes: biphasic media, monophasic 

163 media, and non-homogenizable media. In the biphasic case, the macroscopic description 

16 4 possesses a distinct velocity field for each phase. In the monophasic case both phases 

165 have the same velocity field and we have a single, hybrid, material. In both biphasic and 
we monophasic models, there is only one pressure. The non-homogenizable case is explained 
167 in Appendix D. 

From here on, we assume 



us which implies that at the microscale, the ratio of viscous stresses to pressure in the liquid 

169 phase is 0(e). This includes two biphasic cases, (A4,V) = (0(e 2 ),0(l)) and (A4,V) = 

™ (0(e 3 ), 0(e~ 1 )), and a related monophasic case. We discuss the significance of constraint 

171 (35) in Section 4.1. 

2.4. Main Results 

172 Before proceeding with the expansions, we state our main results. The dependent 

173 variables are V s , the leading order velocity in the matrix, P, the leading order (fluid) 

174 pressure, and , the leading order mean velocity in the fluid. Full notation is given in 




(35) 



175 



Table 5. 
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176 The following systems of equations are derived in Section 3 and Appendices A-C by mul- 

177 tiple scale expansions. They employ an additional assumption that the cell microstructure 

178 possesses certain symmetries which are discussed in Section 3.6 and Appendix E. 

Biphasic-I: In the first biphasic case, V = 0(1) and M. = 0(e 2 ), the leading order 
non-dimensional equations are: 



179 Again, we emphasize that the assumption V = 0(1 = e°) does not imply that the melt and 

wo solid velocities are equal, simply that the ratio of the viscosities is of significantly different 

181 order than the ratio of the velocities. The Im terms are summed over all pairs of I and m. 

182 /c cff is a scalar permeability (and more generally a second order tensor), ( cS . is an effective 

183 scalar bulk viscosity, and for each pair of indices (I, m), rf^ is the anisotropic contribution 

184 to the effective shear viscosity, a second order tensor. These material properties, defined in 
las terms of microscale "cell problems" have been simplified through the domain symmetries. 
186 Derivatives are taken with respect to the dimensionless macroscopic scale x, which we 
is? have suppressed as a subscript for clarity. We note that the equations for the Biphasic-I 
las scaling are in good agreement with previous formulations. The chief difference is the 
189 appearance rj e g, term capturing the grain scale anisotropy, which is new. 



o = p g -vp + v f Ceff. - - 0) J v-v 



(36a) 



+ V • [2(1 - 0)/i s e(V*)] + V • [2 V l ™e lm (y s )} 
0( V '-V) = -^(VP-g') 



(36b) 



V • [<pV f + (1 - 0)V S ] = 



(36c) 
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Biphasic-II: In the second biphasic case, V = 0(e~ l ) and M. = 0(e 3 ), the macroscopic 
system is: 

= pg - VP + V 



Ceff. - 3^(1 - <t>) ) V • V 



(37a) 



+ V • [2(1 - 0)/i s e(V s )] + V • [2rfe m (V s )] 



V > = -^(VP + g') (37b) 
V • (0V 7 ) = (37c) 

wo C e ff., ?7 e ff., and fc e ff. are as above. The first equation is the same as (36a) from the Biphasic- 

191 I model. The differences of the other two equations from (36b - 36c) reflect that when 

192 the fluid velocity is sufficiently greater than the solid velocity and the viscosities are 

193 sufficiently different, the coupling between phases has weakened. In this scaling regime 

194 the two phases only communicate through the pressure gradient. 

Monophasic: In the limit that the melt becomes disconnected, Biphasic-I limits to a 
monophasic system: 

= pg - VP + V • [2p s (1 - 0) e(V s ) + 2?fe m (V s )] (38a) 

V • V s = (38b) 

195 where r) e s, is as above. This is an incompressible Stokes system modeling a composite 
we material with an anisotropic viscosity. 

3. Detailed Expansions and Matching Orders 

197 This section and Appendices A-B provide the detailed derivation and expansions re- 

198 quired to derive the Biphasic-I model summarized in Section 2.4. The other two models 

199 are derived in Appendix C. This material is admittedly technical but we want to provide 

200 sufficient information to offer a road map for related studies. 

DRAFT September 4, 2009, 4:07pm DRAFT 
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We begin by writing our equations in dimensionless form. For the scaling regimes we 
study, the dimensionless forms of the force equations, (11), the in compressibility equations, 
(5b), and the boundary conditions, (6a) and (6b), are: 

V y ■ [-fl + 2e- 1 li s e y (v s )] + eg s = (39a) 

V y ■ [-p f I + 2e 1 ji f e y (v f )} + eg f = (39b) 

V y ■ v s = (39c) 

V y • v f = (39d) 

[-p s I + 2e~ 1 jl s e y (v s )] ■ n = [-p f I + 2e 1 jl f e y (v f )] ■ n (39e) 

v s = Vv / (39f) 

All dependent variables are functions of both x and y. Periodicity in y is imposed to 
capture the periodicity of the microstructure. Derivatives act on both arguments, 



Odd 

dyi dyi dxi 



Analogously, divergence, gradient, and strain rate operators become: 

V„- i-> V y ■ +eV x - (41a) 
Vy^V y + eV, (41b) 
e y I— > e y + ee x (41c) 

3.1. Hierarchy of Equations 

Expanding all variables using Eq. (8) and applying the two scale derivatives, we arrive 
at two hierarchies of equations, one for each phase, which can be solved successively 
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Details of these expansions are given in Appendix A. For the matrix, each iterate is: 

0(e n ): V y ■ a s ^ + V, • a< n ^ + 5 n , lg s = in Y s (42a) 

0(e n+1 ) : V y ■ v s(n+1) + V, • v s(n) = in Y s (42b) 

0{e n ) : (j s(n) • n = a f{n) • n on 7 (42c) 

a s(n) = _ p s{n)j + 2fls [ ex ( v «W) + e y (^ n+1) )] (42d) 

201 where 5 n ^i is the Kronecker delta, such that gravity only acts at order n — 1. Gravity does 

202 not participate in the earlier iterates. Treating a s<n ~ 1 - ) and v s( - n ) as known, the equations 

203 can be interpreted as an inhomogeneous Stokes system for v s ( n+1 ) and p s ( n \ The first 

204 iterate of this system is at n = —1, and we set <r s< ~ 2 ) = v s(_1 ) = cr^ -1 ) = p s ( _1 ) = 0. 

We note that the above equations can be interpreted at each order as a linear system 
in the spirit of the linear algebra problem Ax = b. As with all such problems, there is 
a solvability condition which must be satisfied. For our system, the constraint can be 
interpreted as follows: to be solvable at order n, the integrated surface stress on the solid 
exerted by the fluid, must match the integrated force felt within in the solid, 



/ a*™ ■ ndS = - f (V x ■ a s ^ + 5 n , lg s ) dy. 



(43) 



The enforcement of (43) separates scales and steers us to the macroscopic system. This 
condition can be derived by integrating (42a) over Y s , invoking the divergence theorem 
on the V y ■ a s ^ term, and applying boundary condition (42c). 
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Complementing the equations for the matrix is a hierarchy of equations for the melt 

phase: 

0(e n+1 ) : V y ■ (7 / ( n+1 ) + V, • <r /(ri) + 5 nA g f = in Y) (44a) 
0(e n ) : V y ■ v f ^ + V x ■ v^"" 1 ) = in y> (44b) 

0(en- v/(«) = { vS(n) -7ifV = 0(l), 
1 j - [v^- 1 ) on 7 if V = OCe" 1 ). 1 j 

a /(n) = _ p f(n)j + 2/i/ ^ y /(n-l)) + 

208 As in the solid case, we treat lower order terms, <j^( ra_1 ) and v^ n ~ 2 ) as known, then solve 

209 for pressure and velocity v^ n_1 ). The first iterate of this system is at n — — 1, and 

210 we S et a^" 1 ) = v^ 1 ) = v^ 2 ) = v^ 1 ) = v s (- 2 ) = 0. 

Again, there is a solvability condition. At each order, the flow of the solid at the 
boundary must balance the dilation or compaction of the fluid: 

f ft ^ (~ f v s(n) -ndS ifV = 0(l) , . 

J Yf \- r 7 v s(n_1) -ndS ifV = 0(e- 1 ) V ; 

211 This can be derived by integrating (44b) over Yf, invoking the divergence theorem on the 

212 V y ■ v-^ n ) term, and applying boundary condition (44c). Both solvability conditions (43) 

213 and (45) will be essential for developing macroscopic effective media equations. 
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3.2. Leading Order Equations 

The leading order equations are the same in the three scaling regimes we examine. From 

(42a), (42b), (42c), and (44a), the leading equations are: 

0(e~ l ) : V„ ■ a s{ - t] = in Y s (46a) 

O(e ) : V„ • v s(0) = in Y s (46b) 

O(e ) : V y ■ a /(0) = in Y f (46c) 

0(e- 1 ) : a s{ - t] • n = on 7 (46d) 

These equations can be solved analytically to show that the leading order matrix velocity 
and melt pressure are independent of the fine scale. 

To solve for the leading order solid velocity v s( W, note that the solid stress is a s ^^ = 
2fj Js e y ('v s ^) from (42d) and multiply (46a) by v s (°) and integrate by parts over Y s , 

Jy s 

= I vf ) af l) n J dS-2^ s I \e y (v^) f dy = 
J 1 Jy s 

Applying the boundary condition (46d), we obtain 

/ |e,(v*(°))| 2 rfy = 0. 
Jy s 

which implies that v s (°) is constant in y, 

v a (°)= v i(0) (x) (47) 

v s (°) automatically satisfies (46b). 

Turning to the fluid, the fluid stress is given by (44d) as = —p^I, thus (46c) 
becomes, 

V, • = Vj/ ■ (-p/ (0) /) = -V/°) = 0. 



DRAFT 



September 4, 2009, 4:07pm 



DRAFT 



X - 22 SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

which implies, 

p/(o) =p /(o)( x ). 



(48) 



3.3. Successive Orders in the Solid Phase 

At the next order (n = 0) in (42a-42c), 

0(e°) : V y • a s(0) = in Y s (49a) 

0(e l ) : V, • v s(0) + V y • v s(1) = in Y s (49b) 

O(e ) : <r s(0) • n = a m • n on 7 (49c) 

From (42d), 

^(o) = _^(o)j + 2fIs (^(D) + ex (y-(0))) 

Solvability condition (43) on the stresses is satisfied because = p- f ^ (x) yielding 

jf (-p /(0) /) • n = 0. 

It is helpful to define the pressure difference between solid and fluid as q = p s ^ —p^°\ 
v s(1 ) and q solve 

V y ■ (-ql + 2fi s e y (v s ^)) =0 in Y s (50a) 

V y ■ v s(1) = -V a • v s(0) in Y s (50b) 

(-gJ + 2 f i s ey(^ 1) )) ■ n = (-2fi s e x (v s(0) )) ■ n on 7 (50c) 

This is an inhomogeneous Stokes problem with the forcing terms V x • v s (°) in (50b) and 
2/i s e x (v s (°)) ■ n in (50c); all forcing terms are independent of y. Because the problem 
is linear, we can solve for each component of the forcing independently. The complete 
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solution is the superposition: 

v s(1) = 2e xM (v^) x lm - (V, • v'<°>) ^ (51) 
q = p.(0) _ p m = 2fisex im (v .(0)) ^ _ ^ (Vx . v ,(o )) c (52) 

217 Summation over Im is implied. For each ordered pair (l,m), there is a velocity, x im , 

218 and pressure, n lm , contributed from the corresponding component of the surface stress on 

219 the solid, 2e x j m (\ s<yQS> ). The velocity £ and pressure ( arise from the dilation/compaction 

220 forcing. x' m , ir lm , £, and ( are defined in Table 6 and full statements of the cell problems 

221 are given in Appendix B. These solve the aforementioned auxiliary, or cell, problems, 

222 which are Stokes boundary value problems posed on Y s . 

223 Cell problems may be interpreted as the unit response of the medium to a particular 

224 forcing. For generic three-dimensional cell geometries, the cell problems lack clear analytic 

225 solutions, and one must resort to numerical computation to understand them. In our 

226 second paper, we survey them numerically. 

227 We make two observations on (52). First, it agrees with models that permit the pressures 

228 to be unequal, as in Scott and Stevenson [1984]; Stevenson and Scott [1991]; Bercovici et al. 

229 [2001a]; Bercovici and Ricard [2003]. It also makes clear that the question of whether 

230 there are one or two pressures in macroscopic models of partial melts is entirely semantic. 

231 There are two pressures, but to leading order each can be expressed in terms of the other. 

232 Second, it captures that part of any pressure jump is due to the macroscopic compaction 

233 of the matrix. Such a relation was also discussed in Spiegelman et al. [2007]; Katz et al. 

234 [2007]. 
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3.4. Macroscopic Force Balance in the Matrix 

235 Though we have solved for v s ^\ p s ^ in terms of v s< -°^ and p- f< -°\ we still do not have 

236 a macroscopic equation relating velocity and pressure. To find such an equation we go 

237 to the next order of equations for the matrix and use the solvability condition, (43), to 

238 constrain them. This constraint becomes our macroscopic equation; we do not actually 

239 solve for and p s ^\ 

At the next order of (42a-42c), the equations are: 

Oie 1 ): V x ■ a s{0) + V y ■ + g s = in Y s (53a) 
0(e 2 ) : V, • v s(1) + V y • v s(2) = in Y s (53b) 
Oie 1 ) : a s ^ ■ n = a f ^ ■ n on 7 (53c) 

^(i) i s gi ve n by (42d): 

^(i) = + 2/is [e x (v s ^) + e y (^)] (54) 

According to our force matching solvability condition (43), 

/ a f ^ -ndS=- [ (V x • <x s{0) + g s ) dy (55) 

J-y JY s 

By stress boundary condition (53c), a s ^ • n = • n on 7, so 

f (y x . a m + g ^ dy = _ I a j(i) . ndS 

JY s J-y 

= [ Vya^dy 

Using fluid momentum equation (44a), Vj, • cr^ 1 ^ = — V x ■ a^°"> — in Yf, hence 

I (V, • a<V) dy + [ (V* • a^) dy 

Jy s X JYf ' (57) 

+ (1-0) g s + 0g' = O 



(56) 
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Commuting the integration and divergence operators, 

-V, [( P si0) }s + (p m )f] 

+ 2^ x .[(e x (v s ^)) s + (e y (v s ^)) s ] (58) 
+ pg = 

where angle brackets {-) s denote volume averages over the solid domain Y s (likewise (■)/ 

over Yf). If we substitute (51) and (52) for p s ^ and v s ^\ then 
= pg - V x p m 

- V. {2fi s e xM (v^) {n lm ) s - fi s (OsV x ■ v^} 

(59) 

+ 2^V a • {(1 - 0) e x (v s ^) + 2e xM (v s( V)(e y (x lm ))s} 

-2fi s V x -{(e y (0) s V x -^} 
240 We now have an equation for v s (°) and p^°\ both functions of x. Multiplying this equa- 

2« tion by P s /L restores dimensions. Again, we note that we did not solve (53a-53c). (59) is 

242 merely the equation that must be satisfied for (53a-53c) to satisfy momentum compati- 

243 bility condition (43). 

3.5. Macroscopic Force Balance in the Fluid 

244 We now seek macroscopic equations for the melt. As in the case of the solid, we must 

245 iterate out to the second order correction and use the solvability condition to obtain a 

246 macroscopic equation. 

We first solve for the first correction, obtaining v^°) and p^ x \ and average them. From 
the hierarchy of fluid equations, (44a - 44d), the fluid equations at this order are 

0(e v ) : V x ■ a m + V y ■ + g f = in Yf (60) 
0(e°) : V y ■ v f ^ =0 in Y f (61) 
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with stress 

<T fW = -pfWl + 2 t ji f e y (v f W) 

and boundary conditions 

V /(0) = {VW on 7 if V = 0(1) 

[0 on 7 if V = Oie' 1 ) V ; 

One of the most relevant scaling regimes for magma migration is Biphasic-I with V = 
0(1) and M. = 0(e 2 ), summarized in Section 2.4. We derive it here. The two other 
systems, Biphasic-II and Monophasic, are similar and presented in Appendix C. 

If we substitute the stresses <7^°) and a^ 1 ^ into (60 - 61), we have 

_ Vxp m _ Wyp f(i) + (°> + g' = mY f 

V y • v /(0) = in Y f 

with boundary condition v-^ ) = v s ^°) on 7. Recall that p^°\ v s( -°\ and are interpreted 
as known, inhomogeneous, y independent quantities forcing v^°) and p^ x \ Since it is 
easier to solve a problem with homogeneous boundary conditions, we define w = v^°) — 
v s (°\ simplifying the above equations into 

-V yP m + fi f V 2 y w = V xP m -g f in Y f (63a) 

V y ■ w = in Yf (63b) 
w = on 7 (63c) 

This is the classic homogenization problem of flow in a rigid porous medium and leads to 
Darcy's Law. It is discussed in many of the cited texts on homogenization, particularly 
Hornung [1997]. 



DRAFT 



September 4, 2009, 4:07pm 



DRAFT 



SIMPSON ET AL.: A MULTISCALE MODEL OF PARTIAL MELTS 

The volume compatibility condition (45) is trivially satisfied since w| 7 = 0, 



X - 27 



= / (V„ • w)dy = I w • ndS = 0. 

JYt J-/ 



'Yf 

As in the case of the solid phase, we solve (63a)-(64c) via cell problems, taking advantage 
of the linearity of the problem. We decompose the right hand side forcing terms in (63a) 
into ei, e 2 and e 3 components, solving in each coordinate, then forming the superposition 
of the three to get the solution. Let q\ k* be y periodic functions solving: 



-V y q l + V^k* = -e, inY) 

■*7 



V y ■ k l = in Y t 



k % = on 7 



(64a) 
(64b) 
(64c) 



ei is the unit vector in the i-th direction. These problems thus measure the unit response 
of the fluid to such a forcing. Using the solutions, 



w = --k* (d XlP ^-g{) 
/«(x, y ) = -^ (d XtP W-g{) 



(65) 
(66) 



Averaging over Yt, we get the macroscopic equation for the fluid, 



(v^°)) 



^ (vV (0) - g') 



(67) 



This is Darcy's Law with buoyancy and in a moving frame. (K)f is the permeability 
tensor. K is the matrix, or alternatively the second order tensor, 



and 



K = [k 1 k 2 k 3 ] 



(K) f = \j Yf kMy J Yf k*dy J Yf k 3 rfy 



(6J 



(69) 
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253 While the leading order solid velocity is y-independent, the leading order fluid velocity 

254 remains sensitive to the fine scale. For a macroscopic description, it can only be defined 

255 as an average flux; this is the Darcy velocity of the fluid. 

This is not yet a closed system. Advancing to the next order of (44a - 44c), we have 

0(e 2 ) : V, • <7 /(1) + V y ■ a f W =0 in Y f (70) 
0(e x ) : V. • + V y • v f{1) = in Y f (71) 
0(e 1 ) : v /(1) = v s(1) on 7 (72) 

The solution must satisfy the volume compatibility condition (45), 

V, • yfW = - [ v s (°) • ndS = (73) 

Combining this with (49b), we get 

V,-[(v^) / + (l-^°>]=0 (74) 

256 This is a macroscopic volume compatibility condition. Equations (59), (67), and (74) now 

257 form a closed system. Dimensions may be restored to (67) by multiplying by V* and (74) 

258 by Vf/L; a factor of £ 2 will appear in front of (K)f, as expected. 

3.6. Symmetry Simplifications 

259 The macroscopic equations can be simplified if we assume that the cell geometry is 

260 symmetric with respect to both reflections about the principal axes and rigid rotations. 

261 Though this is a further idealization, the equations retain their essential features. 

Under these two assumptions, (67) (for Biphasic I) and (CI) (for Biphasic II) are 

{vfl% - 0v*(°) = -— (W (0) - g f ) (75) 

H f 

(v^) / = -^(W ( ° ) -g / ) (76) 
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(59) becomes 



= pg - V xP m + V, (Veff. - - 0)) V, • 

+ V, • [2(1 - <j>)p s e x (^) + 2 V l ™e xM (v s ^)} 



(77) 



k e s., Cefr.) an d r/ eS are defined in terms of the solutions of the cell problems: 



fceff. = {K n )f 

2 

Ceff. = fis{C)s ~ oVsO- ~ 0) 



(78) 



(79) 



(80) 



262 ?] cff . is a fourth order tensor. It is a supplementary viscosity, capturing the grain scale 

263 anisotropy of the cell domain. With these symmetry reductions, there are now only four 

264 material parameters to be solved for: (Kn)f, (Q s , (e 2/) n(x 11 )) s , and (e ?/i i2(x 12 ))s corre- 

265 sponding to the macroscopic permeability, bulk viscosity and two effective components of 

266 an anisotropic viscosity. Additional details of the symmetry simplifications may be found 

267 in Appendix E. 

268 If we now define V s = v s (°\ \- f = (v^^)//0, and P = p^°\ and drop the x subscripts 

269 from the derivatives, the above equations become (36a - 36c) presented in Section 2.4. 

4. Discussion 

270 We have successfully used homogenization to derive three macroscopic models for con- 

271 servation of momentum in partially molten systems. We now consider these models fur- 

272 ther, compare them with previous models derived using multiphase flow methods and 

273 discuss some caveats and future directions. 

4.1. Remarks on Homogenization Models 
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The differences amongst the three models of Section 2.4 arise from the assumptions 
on two dimensionless numbers, V and At, and the microstructure. All three rely on the 
additional assumptions that Q{ = 0(e) and V = 0(1). It is helpful to write the three 
models as a unified set of equations: 



274 As V varies from 0(e°) to 0(e _1 ), we transition between Biphasic-I and Biphasic-II. 

275 Letting the pore network disconnect, k e g. — > 0. Consequently, V' — > V -1 V S in (81b). 

276 This recovers macroscopic incompressibility in (81c), V • V s = 0. The divergence terms 

277 also drop from the matrix force balance equation. Making rigorous mathematical sense 

278 of the transition between the connected and disconnected pore network is an important 

279 open problem. It is also interesting that the scalings do not fully describe the macroscopic 

280 equations; the grain scale structure can play a role. 

281 We return to our motivating problem, partially molten rock in the asthenosphere. As 

282 we saw in Section 2.3, for a given e, the parameters V and M. include a range where 

283 a macroscopic description is possible. We lose our ability to homogenize when either 

284 MV 3> e 2 or MV <C e 2 . There may be interesting transitions here. That the two 

285 parameters must be related by MV = 0(e 2 ) would seem a serious constraint on this 

286 approach and its applicability; however, this has another interpretation. 
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(81a) 



+ V • [2(1 - 0)/i,e(V s ) + 2?fe m (V s )] 
0( V /-V- 1 V s ) = -^(VP + g / ) 



(81b) 



V- [0V' + V -1 (l -0)V S ] =0 



(81c) 
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The condition on MV stipulates that the length scales, viscosities, and velocities, be 
related by 

L = V|£ < 82 > 

This also assumes V = 0(1). This can be reinterpreted as the macroscopic length scale 
on which, given the viscosities and characteristic velocities of a partially molten mix we 
should expect to observe a biphasic, viscously deformable, porous media. Based on our 
estimates on the viscosities, velocities, and grain scale in Table 3, 

L « KT 1 - 10 5 km (83) 

Length (82) is similar, but not identical to the compaction length of McKenzie [1984], 



Smm = Mi - m + M 
V 

287 The general scaling is similar as k oc £ 2 therefore the leading scaling is £y/fj^/Jif. Nev- 

288 ertheless <5ms4 is porosity dependent through both permeability, k and the viscosities, ( s 

289 and fx s , making it dynamically and spatially varying. To understand the variation in 

290 compaction length, it is critical to calculate both permeabilities and viscosities that are 

291 consistent with the underlying microstructure. Homogenization provides this computa- 

292 tional machinery through the cell problems. The companion paper calculates consistent 

293 constitutive relations for several simple pore microstructures and suggests that in the limit 

294 — > 0, 5ms4 —* which has important implications for the transition to melt-free regions. 

295 However, L is not a substitute for <5ms4; such a subsidiary length scale may also appear. 

Under the assumption that V = 0(1), (82) also bears resemblance to the compaction 
length of Ricard et al. [2001], 



#brsoi = J^fi (85) 
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296 k,q is a geometric prefactor in a power law scalar permeability relationship k = Ko4> n and 

297 KQ OC £ 2 . 

4.2. Comparison with Existing Models 

298 There are several interesting and important differences between our results and previous 

299 models derived using multiphase flow methods. Most fundamental is that we begin with 

300 a grain scale model, assume certain scalings, and formally derive a macroscopic model. 

301 The anticipated constitutive laws also emerge from these assumptions. 

In the limit of large viscosity variations, the conservation of momentum equations in 
previous models can be closely identified with the Biphasic-I model, where V = 0(1) and 
M. = 0(e 2 ), given by equations (36a - 36c), providing some validation. Compare with 
McKenzie [1984], 

d t (pf(f>) + V • (p/0V^) = mass transfer (86a) 
d t [p s (l - <p)] + V • [p s (l - <f))V s ] = -mass transfer (86b) 
<P(V f - V s ) = - — (VP-gO (86c) 



= pg - VP + V • [2(1 - 0Ke(V s )] + V 



(i-0)(C-^ s )v-v s 



(86d) 



302 /£, and £ s are the permeability, shear viscosity, and bulk viscosity, which have unspeci- 

303 tied dependencies on porosity. We have reused the symbols V^, V s , and P, to denote the 

304 macroscopic fluid and solid velocities, and pressure. 

305 In the absence of melting and freezing, there is good agreement between the two models 

306 if we make the identifications Ceff. = Cs and fc e ff. = The main difference is the appearance 

307 ?7 e ff. term in (36a), reflecting our consideration of a microstructure. We emphasize that 
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08 this macroscopic anisotropy is geometric in origin; the grain scale model was isotropic in 

09 each phase. 

Now we compare with Bercovici and Ricard [2003], in the "geologically relevant limit" 
described by the authors in their Section 3.1. With a bit of algebra, and using our 
notation, this can be written as: 

d t (f) + V • (0V-O = (87a) 
«9 t (l-0) + V-[(l-0)V s ] = O (87b) 
0( V /-V S ) = - — (VP-gO (87c) 
= pg - VP + V • [2(1 - 0)/i s e(V s )] + V 

+ V (surface energy and damage) 

10 In this model, Co is a dimensionless, 0(1) constant. The surface energy and damage 

11 terms, which we have not reproduced, capture surface physics and grain deformation. In 

12 the absence of these physics, there is again good agreement between Biphasic-I and this 
o model if we make the identifications Ccff. = HsCo'fi' 1 and A; cff . = k. As with the McKenzie 
« model, the principal difference comes from the r] e s, term. Bercovici and Ricard [2003] 
is noted that if one eliminates mass transfer in (86a - 86d) and surface physics from (87a - 

16 87d), the two models are identical subject to the identification ( s = /i s Co0 _1 . 

17 The microscale model we homogenized, assuming only fluid dynamical coupling between 
is the phases, was sufficient to generate macroscopic equations consistent with previous 

19 models in the absence of grain-scale surface energies. An important open problem is to 

20 find a grain scale model amenable to homogenization, that includes grain scale diffusion. 

21 One might then see a consistent macroscopic manifestation of these physics, which could 



(l-0)(/^-j^)V-V s 
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322 be compared to models that have already attempted to included them [e.g., Ricard and 

323 Bercovici, 2003; Hier-Majumder et ai, 2006]. 

324 As mentioned previously, an advantage of the homogenization derivation over the mul- 

325 tiphase flow derivation is that there is not the same need for closures. In other models, 

326 one may posit and then seeks closures for permeability, bulk viscosity, shear viscosity, and 

327 interphase force. These parameters might be constrained by other information; however, 
323 this will not yield an inherently self-consistent model. One particularly difficult closure is 

329 the interphase force, the force that one phase exerts on the other. The interphase force, 

330 which is a macroscopic re-expression of the melt-matrix boundary conditions, is poorly 

331 constrained and non-unique. Indeed, the model in Bercovici et al. [2001a], using one inter- 

332 phase force could not replicate the model of McKenzie [1984] . An equally valid interphase 

333 force led to (87a - 87d). As noted, taking out the additional physics, this agrees with 

334 (86a-86d). Though this a desirable result, the non-uniqueness of the terms remains an 

335 issue. 

4.3. Some Caveats 

336 Homogenization provides a more rigorous method for derivation of macroscopic equa- 

337 tions as well as a clear mechanism for computing critical closures. Nevertheless, it is 

338 not foolproof and includes its own set of assumptions whose consequences need to be 

339 understood. 

For example, if the cell domains of Section 2.1 are independent of x, then the porosity 



is constant: 
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340 But a perfectly periodic microstructure is unrealistic. Furthermore, once motion begins, 

341 the interface moves, likely breaking the periodic structure. If the domains do have x 

342 dependence, Yf = i/(x), then we can have = 0(x). This introduces technical difficulties 

343 in (53a), as additional terms for gradients with respect to the domain should now appear. 

344 See Appendix F for details. 

345 A similar omission has been made in the poro-elastic literature [see Lee and Mei, 

346 1997a, b, c; Lee, 2004, for a discussion]. As the elastic matrix deforms, the interface 

347 moves, changing the cell geometry. Earlier work Auriault [1991a]; Hornung [1997]; Mei 

348 and Auriault [1989] implicitly assumed that this deformation was small compared to the 

349 grain scale and could be ignored. This issue also bedevils the sintering and metallurgy 

350 papers Auriault et al. [1992] and Geindreau and Auriault [1999]. In high temperature, tex- 

351 turally equilibrated systems as might be expected in the asthenosphere, grain-boundary 

352 surface forces may help to maintain the geometry of the micro-structure even during 

353 large deformations. However, a consistent homogenization would need to include these 

354 additional microscale processes. 

355 Despite this obstacle, our equations are still of utility in several ways. The first is 
sse that they are a macroscopic description of a constant porosity piece of material. Such a 

357 description has not been rigorously derived before for partially molten rock. It also acts as 

358 a tool for verifying the multiphase flow models. Taking <fi to be instantaneously uniform, 

359 such a model should reduce to our equations. Under simplifications, the other models, 

360 such as McKenzie [1984]; Bercovici and Ricard [2003], are in agreement, up to the r) e Q. 

361 expression. 
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362 Another interpretation is that our models are valid when porosity varies sufficiently 

363 slowly. Under such an assumption, the omitted terms would be higher order in e and 

364 could be justifiably dropped. There is a certain appeal to this; it would not make sense to 

365 discuss the homogenization of a material in which there were tremendous contrasts in the 

366 porosity over short length scales. Moreover, the typical porosity is 0(1%), so that if the 
se? porosity parameter were also scaled, these terms may indeed be small. Such assumptions 

368 of slowly varying porosity underly all of the multiphase flow derivations (and general 

369 continuum mechanics approaches). 

Our final interpretation is that the equations are part of a hierarchical model for partial 
melts. If we ignore melting and assume constant densities, conservation of mass can be 
expressed as 



370 We might then assume that the grain matrix may be approximated by some periodic 

371 structure at each instant. This is consistent with observations. Although the matrix 

372 deforms viscously, it retains a granular structure. Our equations are then treated as the 

373 macroscopic force balances to determine V s , and the system evolves accordingly. 

374 Another issue with the homogenization approach is that though it illuminates how the 

375 effective viscosities and the permeability arise through the cell problems, calculating the 

376 relationship between ( e g,, r) e s, and /c e fr. and the microstructure (e.g. porosity), requires nu- 

377 merically solving the cell problems. The companion paper, Simpson et al. [2008b] , explores 

378 this, calculating effective constitutive relationships for several idealized pore geometries. 
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4.4. Open Problems 

379 There are several ways this work might be extended. A natural continuation is to 

380 model the partial melt as a random medium. This might more realistically model the 

381 pore structure of rocks. The equations for upscaling could also be augmented by giving 

382 the matrix a nonlinear rheology, as in Auriault et al. [1992]; Geindreau and Auriault [1999]. 

383 This may be particularly important for magma migration; a nonlinear matrix rheology 

384 was needed to computationally model physical experiments for shear bands in Katz et al. 

385 [2006] and in general, non-linear power-law rheologies are expected in the dislocation creep 
sse regime [e.g., Hirth and Kohlstedt, 1995b]. 

38? Important mechanisms absent from these equations are surface physics which. In a fluid 

388 dynamical description, these might take the form as surface tension and diffusional terms. 

389 Such terms were posited in the models of Ricard et al. [2001]; Hier-Majumder et al. [2006]; 

390 Bercovici and Ricard [2005, 2003]; Bercovici et al. [2001b, a], but it remains to be shown 

391 how such terms in the macroscopic equations might arise consistently from microscopic 

392 physics that includes grain-scale diffusion and/or mass transfer. 

393 The most serious question remains how to a properly study a medium with macroscopic 

394 and time dependent variations in the structure. This would have implications for the 

395 many physical phenomena that also have evolving micro structures. Recent work in Peter 

396 [2007a, b, 2009] on reaction-diffusion systems in porous media may be applicable. 

Appendix A: Details of the Expansions 

The multiple scale expansions of (8) are applied to p s , pf , v s , and substituted into 
the dimensionless equations (39a - 39f), along with the two scale derivatives, (41a - 41c). 
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Dropping tildes, both melt and matrix strain rate tensors expand as: 

e,(v) ~ e° WvM)] + [e x (v<°>) + e,(v«)] + e 2 [e,(yW) + e y (^)} + . . . 

= e e (0) +e l e (D +e 2 e (2) + _ 

The stress tensors become: 

<x s ^ e- 1 [2/i s e s (°)] + e° [-p'^I + 2/i se s «] + e 1 [-p^I + 2// s e s < 2 >] + . . . 

= e - V(-i) + e cV(o) + e i a -(i) + 
^ ^ e o [_ p /(o) 7 ] + e i [_ p /(i) 7 + 2/i /e ^°)] + e 2 [-/ (2) / + 2/i/e^ 1 )] + . . . 

= a W) + e V/(l) + ^(2) + 

Matching powers of e in equations (A2a) and (A2b) in (39a) and (39b) 
e- 1 V, • a s ^ + e° (V x • + V, • <7 a <°>) + e 1 (V a • a s{0) + V, • + g s ) 

+ ... = 

e V, • + e i (v, ■ ^(°) + V, • a™ + g') + e 2 (V* • a™ + V, • a^) 



(Al) 



(A2a) 



(A2b) 



(A3a) 



(A3b) 



+ ... =0 

Analogously, we substitute the expansions into the incompressibility equations (39c- 
39d), to get 

e°V y ■ v*W + e 1 (V x • v*°> + V, • v^ 1 )) + . . . = 

(A4) 

e°V, • v s (°) + e 1 (V, • v*(°) + V, • v s «) + . . . = 

397 The leading order equations of (A4), V y ■ v-^ ) = and V y ■ v s (°) = 0, reflect that at the 

398 grain scale, both phases are incompressible. 

Making the same power series expansions in the boundary conditions, continuity of 
normal stress, (39e), is 

e" V*- 1 ' ■ n + e° (<j'W - a^) • n + . . . = 0. (A5) 

When V = 0(1), the velocity boundary condition within the cell is 

e o (vS (o)_ v /(o) ) + e i (v ,(i)_ v /(i) ) + > ^ = onT (A6) 
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In this case, the velocities are matched at all orders of e. If instead V = 0(e _1 ), then 

e -l v /(0) + £ ( v /(l) _ v s(0)) 

(A7) 

+ e l( v /(2)_ v -(l)) + ... = onT 

In contrast to the V = 0(1) case, the leading order fluid velocity is independent of the 
solid, and there is cross coupling across orders. 

Appendix B: Cell Problems in the Matrix 

In general, there are two classes of cell problems associated with the matrix phase, and 
a total of seven cell problems. Domain symmetry can reduce the number of unique cell 
problems. 

Bl. Cell Problem for Dilation Stress on Solid 

This addresses the term V T -v s (°) in (50b). This is a less common Stokes problem, with 
a prescribed function in the divergence equation. They are briefly discussed in Temam 
[2001]. Let £, C be y periodic functions solving 

V y ■ (-C/ + 2e„(0) =0 in Y s (Bla) 
V y -l=l inY s (Bib) 
(-(I + 2e y (0) -n = on 7 (Blc) 

The solution measures the response of a unit cell of the matrix to the divergence condition 
(Bib). 

B2. Cell Problem for Surface Stresses on Solid 
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This problem tackles the boundary stress in (50c). Let x im , n lm be y periodic functions 
solving 

V y ■ (-7r /m / + 2e y (x lm )) =0 in Y s (B2a) 
V y • x lm = in Y s (B2b) 
(-7r Zm (% + 2e ytij (x lm )) n j = ~\ Wjm + S im Sji) rij on 7 (B2c) 

we ix lm i 7r ' m ) measure the response of a unit cell of the matrix to a given unit surface stress, 

«7 depending on indices (l,m). Observe that because the tensor on the right hand side of 

408 (B2c), operating on n, is symmetric, the solution to problem (l,m) is the same as the 

409 solution for problem (m,l). For general domains, there are thus six unique cell problems 

410 associated with surface stress forcing. 

Appendix C: Additional Scaling Regimes 

411 In addition to the Biphasic-I regime which we derived in Section 3.5, we presented two 

412 other cases in Section 2.4. These are Biphasic-II, where V = 0(e _1 ) and M. = 0(e 3 ), and 

413 Monophasic, where assumes V = 0(1) and M = 0(e 2 ) and additionally assumes the melt 

414 network is disconnected. Their derivation is given in the next two sections. 

CI. Biphasic-II: Unequal Velocities at the Interface 

In this case V = 0(e _1 ) and M. = 0(e 3 ). Following the scheme of Section 3.5, the 
macroscopic equations are 

(v^) / = -^(W (0) -g / ) (CI) 

V x • = 0. (C2) 

415 Multiplying (CI) by V- f and (C2) by V^/L restores the dimensions of these equations. 
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C2. Monophasic: Magma Bubbles 

416 As in Biphasic-I, we take V = 0(1) and M. = 0(e 2 ). However, we now assume that the 

417 fluid is not topologically connected. The equations are the same at all orders of e as those 

418 appearing in Section 3.5. 

Under this assumption on the microscopic geometry, the permeability cell problems, 
(64a-64c), can be shown to have trivial solutions. k l = for % = 1,2,3, so {k) / = 0. 
Because the melt is trapped it must migrate with the matrix, 

v^°)(x,y) = v s ^(x) (C3) 

Combining (C3) with (74), recovers the incompressibility of the matrix, 

V,, • v s (°) = (C4) 

Dropping the divergence terms from (59) completes the system: 
= pg - V x p m - V x [2fi s e xM (v m )(K lm )s] 

(C5) 

+ 2p s V, • [(1 - 0) ex (v-(°>) + 2e xM (v s(0) )(e y (x lm ))s\ 
4w This is a homogenized incompressible Stokes system for a hybrid material with isolated 

420 very low viscosity inclusions. 

Appendix D: Non-Homogenizable Regimes 

When either Q\ 3> e or Q{ <C e, the system is non-homogenizable. By this we mean that 
it is not possible to upscale equations that faithfully preserve our physical assumptions. 
For instance, if Q\ = 0(1) the pressure gradient balances the viscous forces in the fluid 
and there is no scale separation. Working out the expansions, the leading order velocity 
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and pressure in the fluid solve: 

V, • [-p« > + 2pt f e y (v^)} = inY, (Dl) 
V y • v /(0) = 0, in Y) (D2) 
v /(0) = 0, on 7 (D3) 

The solution is v-^°) = 0. Therefore, 

= yf (v/W + 6 yfW + . . .) (D4) 
= e y/ ( v /(i) + . . .) 

421 This implies that | v^' e | = O(eV^), contradicting our physical assumption that |v^ e | = 

422 0(Vf). While this is mathematically reasonable, the model is unable to produce macro- 

423 scopic fluid velocities of order V*. Other upscaling techniques may succeed here, but 

424 homogenization will not. 

Suppose instead Q{ = 0{t 2 ) or smaller. The fluid equations are then: 

0(e°) : - V y p m = in Yf (D5) 
0(e l ) : - V y p f(1) - V x p m + g f = in Yf (D6) 

The first equation implies pf^ = p-^°)(x). Since V x p^ ^ and are independent of y, 
V y pf^ must also be independent. Since it is periodic in y, it is zero. But this implies 

-V x p m + g / = (D7) 

425 The leading order macroscopic pressure gradient plays no role in balancing the viscous 

426 forces in the solid. This contradicts our assumption that there is always a leading order 

427 non- hydrostatic pressure gradient. 
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Though our assumption on the non-hydrostatic pressure gradient may seem arbitrary, 

there is another important reason to identify cases without such a pressure as non- 
homogenizable. There are problems of interest where gravity plays little role, such as 
Spiegelman [2003]; Katz et al. [2006]. In these cases, would be absent from our equa- 
tions, including (D7). Hence, 

V = ?L V x (pfW + epfW + ...) 



e^-V x (pW + . . .) 



L 

= o(« T ) 



428 This implies the macroscopic fluid pressure gradient is not 0(Pf /L), as hypothesized. 

Appendix E: Cell Problem Symmetries 

429 Let us assume our cell domain is symmetric with respect to the principal axes and 

430 invariant under rigid rotations. This permits simplifications of some of the cell problems. 

In the Darcy cell problem, the off-diagonal entries become zero while the diagonal entries 
are all equal. Thus: 

hs. = (ki)/ (El) 

431 For the surface stress problems, when I ^ to, {^ lm ) s = 0. Only the I, to and to, I entries 

432 of the tensor (e y (x lm )) s are non-zero. For I = to, (n ll ) s = |(1 — 0) and only the diagonal 

433 entries of {e y (x ll ))s are non-zero. The trace of all {e y (x lm ))s tensors is zero. More can be 

434 said about e y (x lm ), but it does not benefit the present analysis. See Simpson [2008] or 

435 Simpson et al. [2008b] for more details. 

436 In the dilation stress problem, the off diagonal terms in (e y (^)) s vanish, and the diagonal 
43? entries are equal to |(1 — </>). 
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Appendix F: Spatial Variation in Cell Domain and Time Dynamics 

If the cells have x dependence, Y} = Yf(x), then it is possible that = 0(x). This 
introduces difficulties in (53a), as terms from gradients with respect to the domain now 
appear. Let us elaborate. For fixed x e fi, we associate a particular cell Y = Y(x), with 
fluid and solid regions defined by the indicator functions 1/ and I s : 

I s :OxFh{0,1} (Fla) 

I j: flxyH{0, 1} (Fib) 

Then 

Y f (x) = {yeY |I / (x,y) = l} (F2a) 

y s (x) = {yey |I s (x,y) = l} (F2b) 



Returning to (53a), 

f V, • a^dy + [ W x ■ a^dy 
Jy 3 Jyj 

= Sy^"' + Sy Vx ' amifdy 

= V X - / a s ^dy- [ a s ^-V x I s dy 

JY, JY 



(F3) 



+ V X - / a f ^dy- [ a f W-V x I f dy 
J Y} Jy 

43a Witness the appearance of the VI S and VI/ terms. This is only an issue for (53a). The 
439 other macroscopic equations remain valid when we allow cell variation. 
A second problem is manifest when we consider time dynamics. 
dt<j> = dt I ldy=[v f - ndS 

JYf JT 

= - j \ s -ndS = - J (v s(0) + ev s(1) + ...)• ndS 
Since v s (°) is independent of y, the first term drops. Substituting (51), 

d t <P = -e Ny v s «dy + 0(e 2 ) = eV x • v s ^\l - 0) + 0(e 2 ) 
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440 To leading order, the matrix can only vary by dilation and compaction. 
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Table 1. Notation for domains in homogenization model. 

Symbol Meaning 

7 Interface between melt and matrix within cell Y 

T Total macroscopic interface between melt and matrix 

Q Total macroscopic space occupied by both melt and matrix 

Q f Portion of macroscopic space occupied by melt 

Q s Portion macroscopic space occupied by matrix 

Y The unit cell 

Yf Portion of unit cell occupied by melt 

Y s Portion of unit cell occupied by matrix 
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Figure 1. The macroscopic domain Q. The matrix occupies the gray region while the 



melt occupies the white inclusions. 



Y 




Figure 2. The cell domain, Y, divided into fluid and solid regions, Yf and Y s . The two 
phases meet on interface T. 
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Figure 3. SEM images of synthetic quartzites and marbles from Figure 5 of Wark and 
Watson [1998]. Similar microstructures are seen in olivine basalt aggregates [e.g. Hirth 
and Kohlstedt, 1995a, b] 
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Table 2. Notation for fields in homogenization model. 



Symbol Meaning 



e(v) 


Strain rate tensor, e(v) = ^(Vv + (Vv) T ) 





Volume fraction of melt, = J Y dy 


g 




g / 


P/g 


g S 


Psg 


pf 


Melt pressure 


pfti) 


Melt pressure at order j in the series expansion 


p s 


Melt pressure 


ps(j) 


Matrix pressure at order j in the series expansion 


a* 


Melt Stress Tensor 


a f(j) 


Melt Stress Tensor at order j in the series expansion 


a s 


Matrix Stress Tensor 


a s{j) 


Matrix Stress Tensor at order j in the series expansion 




Melt velocity 


v f(j) 


Melt velocity at order j in the series expansion 


V s 


Matrix Velocity 




Matrix velocity at order j in the series expansion 
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Table 3. Notation and measurements for models of partial melts. 



Symbol 


Meaning 


Value 




Volume Fraction of Melt 


.01%- 10% 


9 


Gravity 


9.8 m/s 2 


£ 


Grain Length Scale 


1 -10 mm 


L 


Macroscopic Length Scale 


1 m - 10 km 


/'/ 


Melt Shear Viscosity 


1-10 Pa s 




Matrix Shear Viscosity 


10 i5_ 10 2i Pa s 


Re/ 


Reynolds Number of Melt 


io~ 8 -io- 5 


Re| 


Reynolds Number of Matrix 


10 -30_ 10 -22 


Pf 


Melt Density 


2800 kg/m 3 


Ps 


Matrix Density 


3300 kg/m 3 


v' 


Characteristic Melt Velocity 


1-10 m/yr 


V s 


Characteristic Matrix Velocity 


1-10 cm/yr 




Yf Y s 

Figure 4. A cell geometry in which both the fluid region, lj, and the solid region, Y s 
are topologically connected. 
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Table 4. Dimensionless numbers for homogenization model. 



Symbol 


Meaning 


Estimate 


e 


Length scale ratio, e = £/L 


o(io- 7 - 10~ 2 ) 


M 


Viscosity ratio, M = /if/pis 


0(l(r 21 - 10" 14 ) 


V 


Pressure ratio, V = P? / P s 


O(e ) 


Ql 


Ratio of viscous force to pressure gradient in 
melt, Q{ = fa f V')/(Pf£) 


0(e 9 - e°) 


e? 


Ratio of viscous force to pressure gradient in 
matrix, Qf = (fi s V s )/{P s £) 


o( e - x ) 




Ratio of viscous force to pressure gradient in 
matrix, Q S L = (^ S V S )/(P S L) 


O(e ) 


n{ 


Ratio of buoyancy force to pressure gradient 
in melt, Tl{ = {p f g£)/P f 




n\ 


Ratio of buoyancy force to pressure gradient 
in matrix, TZj = (p s q£)/P s 


Oie 1 ) 




Buoyancy force to pressure gradient ratio in 
matrix, 1Z S L = (p s gL)/P s 


O(e ) 


V 


Velocity ratio, V = V f /V s 


o(io- 7 - io- 2 ) 
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Table 5. Effective quantities derived by homogenization. 
Symbol Meaning 



77 e ff. Effective supplementary anisotropic viscosity, a fourth order tensor. 
r/ l e Q = 2fi s (x lrn ) s is a second order tensor. 

(K) f Permeability, a second order tensor. The i-th column is given by 

k e fi. Isotropic permeability. Under symmetry, k e e, = f 

(■)/ Volume average of a quantity over the melt portion of a cell, (•)/ = 

J Yf - d y 

(■)s Volume average of a quantity over the matrix portion of a cell, 

(•)» = I Ys - d y 

P Effective macroscopic (fluid) pressure, P = p^°\ 
V / Effective macroscopic fluid velocity, = (v-^°))//</>. 
V s Effective macroscopic solid velocity, V s = v s (°). 
Cefr. Effective bulk viscosity of the matrix, £ e fj. = fi s (C)s — |/Us(l — <P) 
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Table 6. Notation for cell problems. 



Symbol Meaning 

X lm Velocity of the cell problem for a unit shear stress forcing on the solid in the 
Im component of the stress tensor 

Unit vector in the i-th coordinate, ej = (1, 0, 0) 

k* Velocity of the cell problem for a unit forcing on the fluid in the direction 

£ Velocity of the cell problem for a unit forcing on the divergence equation 

n lm Pressure of the cell problem for a unit shear stress forcing on the solid in the 
Im component of the the stress tensor 

qi Pressure of the cell problem for a unit forcing on the fluid in the ej direction 

( Pressure of the cell problem for a unit forcing on the divergence equation 



Table 7. Additional notation for the other models. 

Symbol Meaning 

C An 0(1) Constant 

k Permeability 

kq Permeability constant for a power law permeability, k = /to0 

( s Bulk viscosity 
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